home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / filters / filter.f < prev    next >
Text File  |  1988-10-09  |  23KB  |  315 lines
  1.       implicit real*8 (a-h,o-z)                                         
  2.       real*4 x1,x2,x3,x4,x5                                             
  3.       common/b/cn(30),cd(30),mn,md,const                                
  4.       call reset                                                        
  5.       write(6,1000)                                                     
  6. 1000  format( ' '/                                                      
  7.      &' f3=0-> low or highpass. f1=passband cutoff. f2=stopband cutoff'/
  8.      &' f1<f2 -> lowpass.'/                                             
  9.      &' f3>0 -> bandpass. f1,f2 are limits of passband. f3 is limit of'/
  10.      &' either high or low stopband._ we require f1<f2.'/               
  11.      &' ripple=passband ripple in db. atten=stopband attenuation in db.'
  12.      &)                                                                 
  13.       write(6,910)                                                      
  14. 910   format( ' the filter coefficients will be found in a file called
  15.      $ fort.7'/
  16.      $' the listing of the program output will be found in fort.8'/
  17.      $ ' enter the following arguments in floating point form.')        
  18.       write(6,900)                                                      
  19. 900   format( ' enter sampling rate.')                                  
  20.       read(5,901) sr                                                    
  21. 901   format(g10.0)                                                     
  22.       xnyq=sr/2.d0                                                      
  23.       write(6,903)                                                      
  24. 903   format( ' enter f1')                                              
  25.       read( 5,901) x1                                                   
  26.       write(6,905)                                                      
  27. 905   format( ' enter f2')                                              
  28.       read( 5,901) x2                                                   
  29.       write(6,906)                                                      
  30. 906   format( ' enter f3')                                              
  31.       read( 5,901) x3                                                   
  32.       write(6,907)                                                      
  33. 907   format( ' enter ripple')                                          
  34.       read( 5,901) x4
  35.       write(6,908)                                                      
  36. 908   format( ' enter attenuation in db')                               
  37.       read( 5,901) x5                                                   
  38. c      write(8,98) x1,x2,x3,x4,x5                                       
  39.       write(6,97) x1,x2,x3,x4,x5   
  40. c98    format( ' f1=',f8.3,' f2=',f8.3,' f3=',f8.3,' ripple=',          
  41. c     $f8.3,' atten=',f8.3)                                             
  42. 97    format('c     f1=',f8.3,' f2=',f8.3,' f3=',f8.3,' ripple=', 
  43.      $f8.3,' atten=',f8.3)
  44.       f1=x1                                                             
  45.       f2=x2                                                             
  46.       f3=x3                                                             
  47.       ripple=x4                                                         
  48.       atten=x5                                                          
  49.       call ellips(f1,f2,f3,ripple,atten,sr)                             
  50.       call fresp(200,sr,0.d0,xnyq,f1)                                   
  51.       m2=mn/2                                                           
  52.       write(7,1234) m2,x1,x2,x3,x4,x5                                   
  53. 1234  format('/*elliptical filter with',i3,' sections.'/                  
  54.      &'f1,f2,f3=',3g10.4,' ripple=',g10.4,' db=',g10.4,'*/')                 
  55.       write(7,899) m2,(cn(i),cd(i),i=1,mn)                                 
  56. 899   format('ell(start,dur,inskip,chin,chout,',i3,',',4(g15.8,','))
  57.       write(7,898) const
  58. 898   format(e15.8,')')
  59.       stop                                                              
  60.       end                                                               
  61.       subroutine reset                                                  
  62.       implicit real*8 (a-h,o-z)                                         
  63.       common/b/cn(30),cd(30),mn,md,const                                
  64.       mn=0                                                              
  65.       md=0                                                              
  66.       write(8,200)                                                      
  67.   200 format('//')                                                      
  68.       do 100 m=1,30                                                     
  69.       cn(m)=0.                                                          
  70.       cd(m)=0.                                                          
  71. 100   continue                                                          
  72.       return                                                            
  73.       end                                                               
  74.       subroutine ellips(f1,f2,f3,ripple,atten,samr)                     
  75. c   designs an elliptic filter. all parameters real*8 .                 
  76. c   f3=0 -> lowpass or highpass. f1=passband cutoff. f2=stopband cutoff.
  77. c   f1<f2 -> lowpass.                                                   
  78. c   f3>0 -> bandpass. f1,f2 are limits of passband. f3 is limit of      
  79. c   either high or low stopband. we require f1<f2.                      
  80. c   ripple=passband ripple in db. atten=stopband attenuation in db.     
  81. c   samr=sampling rate in hz.                                           
  82. c    after gold+rader; written by bilofsky, revised by steiglitz        
  83. c    pp.61-65 (elliptic filters), 72,76 (mappings                       
  84. c    from s-plane to z-plane), 87 (approximation                        
  85. c    for u0 and evaluation of elliptic functions).                      
  86.       implicit real*8 (a-h,o-z)                                         
  87.       real*8 k,k1,kay,kprime,k1prim ,nn,kk,kkp,kk1,kk1p                 
  88.       common/ellipt/k,kprime,cosp0,w1,hpass                             
  89.       prime(dummy)=dsqrt(1.d0-dummy**2)                                 
  90.       bpt(w)=dabs((cosp0-dcos(w))/dsin(w))                              
  91.       pi=3.14159265358979d0                                             
  92.       w1=2.d0*pi*f1/samr                                                
  93.       w2=2.d0*pi*f2/samr                                                
  94.       w3=2.d0*pi*f3/samr                                                
  95.       hpass=0.d0                                                        
  96.       cosp0=0.d0                                                        
  97.       if(f3.gt.0.d0)goto1                                               
  98.       if(f1.lt.f2)goto2                                                 
  99. c  modify frequencies for high pass.                                    
  100.       w1=pi-w1                                                          
  101.       w2=pi-w2                                                          
  102.       hpass=1.d0                                                        
  103. c  compute analog frequencies for low/high pass                         
  104.     2 w1=dtan(.5d0*w1)                                                  
  105.       w2=dtan(.5d0*w2)                                                  
  106.       goto3                                                             
  107. c  compute analog frequencies for band pass.                            
  108.     1 cosp0=dcos((w1+w2)/2.d0)/dcos((w1-w2)/2.d0)                       
  109.       w1=bpt(w1)                                                        
  110.       de=w3-w2                                                          
  111.       if (de.lt.0.d0) de=w1-w3                                          
  112.       w2=dmin1(bpt(w1-de),bpt(w2+de))                                   
  113. c  compute params for poles,zeros in lambda plane                       
  114. 3     k=w1/w2                                                           
  115.       kprime=prime(k)                                                   
  116.       eps=dsqrt(10.d0**(.1d0*ripple)-1.d0)                              
  117.       a=10.d0**(.05d0*atten)                                            
  118.       k1=eps/dsqrt(a*a-1.d0)                                            
  119.       k1prim =prime(k1)                                                 
  120.       kk=kay(k)                                                         
  121.       kk1=kay(k1)                                                       
  122.       kkp=kay(kprime)                                                   
  123.       kk1p=kay(k1prim )                                                 
  124.       n=idint(kk1p*kk/(kk1*kkp))+1                                      
  125.       nn=n                                                              
  126.     5 u0=-kkp*dlog((1.d0+dsqrt(1.d0+eps*eps))/eps)/kk1p                 
  127. c  now compute poles,zeros in lambda plane,                             
  128. c    transform one by one to z plane.                                   
  129.       dd=kk/nn                                                          
  130.       tt=kk-dd                                                          
  131.       dd=dd+dd                                                          
  132.       n2=(n+1)/2                                                        
  133.       do 4 i=1,n2                                                       
  134.       if (i*2.gt.n) tt=0.d0                                             
  135.       call stuff1(-kkp,tt,'zero')                                       
  136.       call stuff1(u0,tt,'pole')                                         
  137. 4     tt=tt-dd                                                          
  138.       return                                                            
  139.       end                                                               
  140.       subroutine stuff1(q,r,whatsi )                                    
  141. c    transforms poles and zeros to z-plane; stuffs coeff. array         
  142.       implicit real*8 (a-h,o-z)                                         
  143.       real*8 k,kprime                                                   
  144.       common/b/cn(30),cd(30),mn,md,const                                
  145.       character*4 whatsi                                                
  146.       complex*16 dcmplx,cdsqrt,dconjg,z,s                               
  147.       common/ellipt/k,kprime,cosp0,w1,hpass                             
  148.       call djelf(snr,cnr,dnr,r,kprime*kprime)                           
  149.       call djelf(snqp,cnqp,dnqp,q,k*k)                                  
  150.       omega=1-snqp*snqp*dnr*dnr                                         
  151.       if ( omega .eq. 0.d0 ) omega=1.d-30                               
  152.       sigma=w1*snqp*cnqp*cnr*dnr/omega                                  
  153.       omega=w1*snr*dnqp/omega                                           
  154.       s=dcmplx(sigma,omega)                                             
  155.       j=1                                                               
  156.       if (cosp0.eq.0.d0) goto 1                                         
  157.       j=-1                                                              
  158.     4 z=(-cosp0+dfloat(j)*cdsqrt(cosp0*cosp0+s*s-1.d0))/(s-1.d0)        
  159.       go to 3                                                           
  160.     1 z=(1.d0+s)/(1.d0-s)                                               
  161.       if(hpass.ne.0.d0)z=-z                                             
  162.     3 if(dabs(dimag(z)).le.10.d-10) goto 2                              
  163.       if(dimag(z).lt.0.d0) z=dconjg(z)                                  
  164.       if(whatsi.eq.'pole')goto5                                         
  165.       mn=mn+1                                                           
  166.       cn(mn)=-2.d0*dreal(z)                                             
  167.       mn=mn+1                                                           
  168.       cn(mn)=dreal(z)**2+dimag(z)**2                                    
  169.       goto6                                                             
  170.     5 md=md+1                                                           
  171.       cd(md)=-2.d0*dreal(z)                                             
  172.       md=md+1                                                           
  173.       cd(md)=dreal(z)**2+dimag(z)**2                                    
  174.     6 write(8,202)whatsi,z                                              
  175.   202 format(' complex ',a4,' pair at ',d17.9,' +-j',d17.9)             
  176.       if(j.gt.0.or.r.eq.0.d0)return                                     
  177.       j=1                                                               
  178.       go to 4                                                           
  179.     2 x=dreal(z)                                                        
  180.       if(whatsi.eq.'pole')goto7                                         
  181.       mn=mn+1                                                           
  182.       cn(mn)=-x                                                         
  183.       mn=mn+1                                                           
  184.       cn(mn)=0.d0                                                       
  185.       goto8                                                             
  186.     7 md=md+1                                                           
  187.       cd(md)=-x                                                         
  188.       md=md+1                                                           
  189.       cd(md)=0.d0                                                       
  190.     8 write(8,201)whatsi,x                                              
  191.   201 format(' real ',a4,' at ',d17.9)                                  
  192.       if(j.gt.0) return                                                 
  193.       j=1                                                               
  194.       go to 4                                                           
  195.       end                                                               
  196.       subroutine fresp(k,samr,f1,f2,f3)                                 
  197. c    plots k pts. of freq. resp. from f1 to f2, norm. at f3             
  198.       implicit real*8 (a-h,o-z)                                         
  199.       complex*16 dcmplx,cdexp,tf,zm,zm2                                 
  200.       common/b/cn(30),cd(30),mn,md,const                                
  201.       pi=3.14159265358979d0                                             
  202.       m2=mn/2                                                           
  203.       write(8,200)m2,(cn(i),cd(i),i=1,mn)                               
  204.   200 format('elliptic filter with ',i5,' sections'/4(d17.9))           
  205.       w=pi*f3/(.5d0*samr)                                               
  206.       zm=cdexp(dcmplx(0.d0,-1.d0*w))                                    
  207.       zm2=zm*zm                                                         
  208.       tf=(1.d0,0.d0)                                                    
  209.       do 1 i=1,mn,2                                                     
  210.     1 tf=tf*(1.d0+cn(i)*zm+cn(i+1)*zm2)/(1.d0+cd(i)*zm+cd(i+1)*zm2)     
  211.       const=1.d0/cdabs(tf)                                              
  212.       write(8,201)const                                                 
  213.   201 format(' const=',d17.9)                                           
  214.       write(8,205)                                                      
  215.   205 format('/   freq     phase',10x,'    amp',10x,'    db.')          
  216.       do 3 j=1,k                                                        
  217.       freq=f1+(f2-f1)*dfloat(j-1)/dfloat(k-1)                           
  218.       w=pi*freq/(.5d0*samr)                                             
  219.       zm=cdexp(dcmplx(0.d0,-1.d0*w))                                    
  220.       zm2=zm*zm                                                         
  221.       tf=dcmplx(const,0.d0)                                             
  222.       do 2 i=1,mn,2                                                     
  223.     2 tf=tf*(1.d0+cn(i)*zm+cn(i+1)*zm2)/(1.d0+cd(i)*zm+cd(i+1)*zm2)     
  224.       amp=cdabs(tf)                                                     
  225.       if(amp.le.1.d-20)amp=1.d-20                                       
  226.       x=dreal(tf)                                                       
  227.       y=dimag(tf)                                                       
  228.       phase=0.d0                                                        
  229.       if(x.eq.0.d0 .and. y.eq.0.d0)goto4                                
  230.       phase=(180.d0/pi)*datan2(y,x)                                     
  231.     4 db=20.d0*dlog10(dmax1(amp,1.d-40))                                
  232.     3 write(8,202)freq,phase,amp,db                                     
  233.   202 format(' ',f10.2,2d17.9,f12.4)                                    
  234.       return                                                            
  235.       end                                                               
  236.       double precision function kay(k)                                  
  237. c    computes kay(k)=inverse sn(1)                                      
  238. c    hastings, approx. for dig. comp., p. 172                           
  239.       implicit real*8 (a-h,o-z)                                         
  240.       double precision k,eta,peta,kk                                    
  241.       dimension a(5),b(5)                                               
  242.       data a/1.38629436112d0, .09666344259d0, .03590092383d0,           
  243.      1    .03742563713d0, .01451196212d0/                               
  244.       data b/.5d0, .12498593597d0, .06880248576d0, .03328355346d0,      
  245.      1 .00441787012d0/                                                  
  246.       kay=a(1)                                                          
  247.       kk=b(1)                                                           
  248.       eta=1.d0-k*k                                                      
  249.       peta=eta                                                          
  250.       do 1 i=2,5                                                        
  251.       kay=kay+a(i)*peta                                                 
  252.       kk=kk+b(i)*peta                                                   
  253. 1     peta=peta*eta                                                     
  254.       kay=kay-kk*dlog(eta)                                              
  255.       return                                                            
  256.       end                                                               
  257.       subroutine djelf(sn, cn, dn, x, sck)                              
  258. c     ssp program: finds jacobian elliptic functions sn,cn,dn.          
  259.       implicit real*8 (a-h,o-z)                                         
  260.       dimension ari(12),geo(12)                                         
  261.       double precision sn,cn,dn,x,sck,ari,geo,cm,y                      
  262.      c,a,b,c,d                                                          
  263.       cm=sck                                                            
  264.       y=x                                                               
  265.       if(sck)3,1,4                                                      
  266.     1 d=dexp(x)                                                         
  267.       a=1.d0/d                                                          
  268.       b=a+d                                                             
  269.       cn=2.d0/b                                                         
  270.       dn=cn                                                             
  271.       a=(d-a)/2.d0                                                      
  272.       sn=a*cn                                                           
  273.     2 return                                                            
  274.     3 d=1.d0-sck                                                        
  275.       cm=-sck/d                                                         
  276.       d=dsqrt(d)                                                        
  277.       y=d*x                                                             
  278.     4 a=1.d0                                                            
  279.       dn=1.d0                                                           
  280.       do 6 i=1,12                                                       
  281.       l=i                                                               
  282.       ari(i)=a                                                          
  283.       cm=dsqrt(cm)                                                      
  284.       geo(i)=cm                                                         
  285.       c=(a+cm)*.5d0                                                     
  286.       if(dabs(a-cm)-1.d-9*a)7,7,5                                       
  287.     5 cm=a*cm                                                           
  288.     6 a=c                                                               
  289.     7 y=c*y                                                             
  290.       sn=dsin(y)                                                        
  291.       cn=dcos(y)                                                        
  292.       if(sn)8,13,8                                                      
  293.     8 a=cn/sn                                                           
  294.       c=a*c                                                             
  295.       do 9 i=1,l                                                        
  296.       k=l-i+1                                                           
  297.       b=ari(k)                                                          
  298.       a=c*a                                                             
  299.       c=dn*c                                                            
  300.       dn=(geo(k)+a)/(b+a)                                               
  301.     9 a=c/b                                                             
  302.       a=1.d0/dsqrt(c*c+1.d0)                                            
  303.       if(sn)10,11,11                                                    
  304.    10 sn=-a                                                             
  305.       goto 12                                                           
  306.    11 sn=a                                                              
  307.    12 cn=c*sn                                                           
  308.    13 if(sck)14,2,2                                                     
  309.    14 a=dn                                                              
  310.       dn=cn                                                             
  311.       cn=a                                                              
  312.       sn=sn/d                                                           
  313.       return                                                            
  314.       end                                                               
  315.